Interplay between phase ordering and roughening on growing films 
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We study the interplay between surface roughening and phase separation during the growth of 
binary films. RenormaUzation group calculations are performed on a pair of equations coupling the 
interface height and order parameter fluctuations. We find a larger roughness exponent at the critical 
point of the order parameter compared to the disordered phase, and an increase in the upper critical 
dimension for the surface roughening transition from two to four. Numerical simulations performed 
on a solid-on-solid model with two types of deposited particles corroborate some of these findings. 
However, for a range of parameters not accessible to perturbative analysis, we find non-universal 
behavior with a continuously varying dynamic exponent. 
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I. INTRODUCTION 

Thin solid films are grown for a variety of technolog- 
ical applications, using molecular beam epitaxy (MBE) 
or vapor deposition. In order to create materials with 
specific electronic, optical, or mechanical properties, of- 
ten more than one type of particle is deposited. When 
the particle mobility in the bulk is small, surface configu- 
rations become frozen in the bulk, leading to anisotropic 
structures that refiect the growth history, and are dif- 
ferent from bulk equilibrium phases. If, for instance, a 
combination of particles are deposited that tend to phase 
separate at the surface, the grown films have lamellae or 
columns of the twophases that extend parallel to the 
growth direction P, Q . This process of phase separation, 
as well as other ordering phenomena, can be affected by 
elastic forces, by the orientation of the growing crystal, 
by properties of the substrate, and by surface roughness. 
The range of possible scenarios is very rich and far from 
understood. There are a variety of analytical and com- 
puter models which try to shed light on some of these 
phenomena, but a systematic study and understanding 
of the possible phase transitions does not yet exist. 

In this paper, we focus on the interplay between phase 
separation and surface roughening, neglecting the possi- 
ble influence of elastic forces, substrate properties, and 
orientation dependencies due to the crystal structure. 
There exist several theoretical studies of phase separa- 
tion during growth that neglect also the effect of surface 
roughness. In all these models it is assumed that the 
mobility of the atoms in the bulk is zero, such that all 
of the dynamics occurs at the surface. A model in which 
the probability that an incoming atom sticks to a given 
surface site depends on the state of the neighboring sites 
in the layer below Q , leads to a phase separation transi- 
tion in the universality class of ordinary Ising models, if 
the model is symmetric with respect to the two phases. 
The same conclusion applies a model in which the top 
layer is fully thermally equilibrated before the next layer 



is added A model for spinodal decomposition during 
growth was introduced in Ref. In this model, phase 
separation is due to surface diffusion, and is limited due 
to the current of incoming particles, leading to a charac- 
teristic scale for the thickness of lamellae or columns, as 
confirmed by Monte-Carlo simulations 6]. 

However, the layer by layer growth mode underlying 
these models is unstable, and growing surfaces generally 
are rough. Several studies exist that investigate growth 
models that contain both phase separation and surface 
roughness. Simulations of an Eden model with two types 
of particles suggest that the surface roughness increases 
due to the phase separation 7]. A solid-on-solid growth 
model where the adsorption probabilities for the two 
types of particles depend on the local neighborhood in 
the layer below leads also to an increased surface rough- 
ness The reason is that particles are more likely to 
be adsorbed within domains than at domain boundaries. 
On length scales much larger than the domain size, a 
crossover to the scaling behavior of the Kardar-Parisi- 
Zhang (KPZ) equation is found. Another computer 
model where particles are adsorbed randomly and subse- 
quently diffuse along the surface leads to domains whose 
thickness is a nonmonotonous function of the deposition 
rate and the temperature, and for a certain range of pa- 
rameter values, the height profile has steep steps at do- 
main boundaries [ic| . A set of coupled Langvin equations 
for this model is suggested in Ref. J,l] and studied using 
stability analysis and Fourier decomposition. 

These studies are rather incomplete, and in particular 
lacking a discussion of the possible effects of the height 
profile on the phase separation dynamics. A first attempt 
to a systematic study of the possible phases and scaling 
behaviors of coupled phase separation and roughening 
during growth was presented in a recent letter by us |l2j| . 
A set of two coupled Langevin equations was suggested, 
and computer simulations in 1+1 dimensions were per- 
formed, revealing a rich phase diagram. It is the purpose 
of this paper to extend and deepen this short study by 
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presenting a renormalization group (RG) analysis, and 
further simulation results. While the RG analysis gives 
information about the behavior of the system in dimen- 
sions close to 4+1 and higher, computer simulations are 
particularly efficient in low dimensions. The general re- 
sults obtained by the two approaches are compatible with 
each other. In addition, the computer simulations in 1-1-1 
dimension reveal interesting nonuniversal behavior for a 
range of parameters that cannot be studied using pertur- 
bative RG. 

The outline of the paper is as follows: In section|nl we 
introduce and discuss the coupled set of Langevin equa- 
tions used in this paper. Scaling laws and critical expo- 
nents will also be defined. Section ITTll presents results of 
the RG analysis of these equations. One of the main find- 
ings is that the lower critical dimension for the surface 
roughening transition is increased from 2 to 4 dimen- 
sions due to the coupling to the critical phase ordering 
dynamics. Section Hvl presents results of computer simu- 
lations. Section analyses the connection of our model 
with the advection of a passive scalar in a velocity field, 
and with directed polymers drifting through a random 
medium. Section IVII contains a summary and discussion 
of our results. 



II. EQUATIONS OF MOTION AND SCALING 
LAWS 

We consider the growth of a binary alloy on a d- 
dimensional substrate. Let x be the coordinate perpen- 
dicular to the growth direction, and t the time. Since 
we assume that all dynamics occurs at the surface of 
the growing material, the equations of motion can be ex- 
pressed in terms of x and t alone. In order to characterize 
surface roughness and phase ordering, we introduce the 
height variable h{x.,t), which is the surface profile at po- 
sition X at time t, and an order parameter m(x, t), which 
is the difference in the densities of the two particle types 
at the surface at position x and time t. The interplay be- 
tween the fluctuations in to, and the height h is captured 
phenomenologically by the coupled Langevin equations, 

dth = :yV'h+^iVhf + ^m^ + a, (1) 
dtm = K{\7^m — rm — um^) + dVh ■ Vto + bm\7^h 
+ ^ni{Vhf + (2) 

with 

(a(x,t)a(x',i')) - 2Dh5\^-^)5(t-t'), 
(C™(x,t)C™(x',i')> - 2A„<5'(x-x')^(i-0- 

Since we are interested in the critical behavior of 
the model, we have assumed that it has the symmetry 
TO — !■ —TO, and included all potentially relevant terms 
compatible with this symmetry. In experiments or com- 
puter simulations, this symmetry can be achieved by tun- 



ing the ratio between the two types of adsorbed parti- 
cles to the appropriate value. In the absence of such an 
order parameter symmetry, the system may undergo a 
flrst-order phase transition which is not considered here. 
Equation is the Kardar-Parisi-Zhang (KPZ) equation 
for surface growth, plus the leading coupling to the 
order parameter. Equation ^ is the time dependent 
Landau-Ginzburg equation for a (non-conserved) Ising 
model, with three different couplings to the height fluc- 
tuations. The Gaussian, delta-correlated noise terms, C/i 
and Cm , mimic the effects of faster degrees of freedom. 

These equations apply to growth by vapor deposition, 
with particles sticking at surface sites with a probability 
that depends on the local environment in the growing 
fllm. The coupling terms in the Langevin equations 
and Q have obvious meaning in this context: The term 
proportional to a implies that particles are more likely to 
be absorbed within domains where they feel a stronger 
attractive force (if a > 0). A negative a can also be 
meaningful: if the adsorption rate within domains is lim- 
ited by the availability of particles of the correct type, 
this can slow down growth. However, if this is due to the 
vapor phase not being well stirred, additional equations 
for the particle concentrations in the vapor phase will be 
needed. Such equations are not included in this paper. 
The contribution from a (with a > 0) implies that do- 
main walls tend to be driven downhill; e.g. if the identity 
of a newly adsorbed particle is more likely to be affected 
by its uphill neighbors than by the downhill ones. A 
positive h indicates that new domains are more likely to 
be formed in hilltops where there are less neighbors that 
could influence the type of particle to be adsorbed. The 
term proportional to c is similar in character to the KPZ 
nonlinearity A, and means that susceptibility to phase 
separation depends on the slope. 

Models for MBE typically assume that particle depo- 
sition at the surface is random, and that no desorption of 
particles takes place. In this case, the height proflle and 
the order parameter dynamics are shaped by diffusion of 
particles along the surface. This physical situation leads 
to a different set of Langevin equations, 

dth = uV'^h+^V'^m^ + Ch, (3) 

dtm = ifV^TO - WTO + Cm, (4) 

where we have again imposed the symmetry to —to. 
The noise terms Cm and Cft, have the same nonconserved 
correlations as for the vapor-deposition model above, due 
to the incoming particle current. Because of the conser- 
vation of volume during surface diffusion, the determinis- 
tic terms on the right-hand side of Eq. ^ must be the di- 
vergence of a current, disallowing the terms proportional 
to A and a in Eq. (^). A negative value of (3 means that 
particles are more likely to be adsorbed within domains 
if they are not needed to the same extent in the neighbor- 
hood. The equation for the order parameter has also the 
form of the divergence of a current (we have only included 
the lowest-order term), plus a nonconserved contribution 
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—vm due to the incoming current of particles that tends 
to reduce the value of the order parameter. The lowest 
order coupling to the height variable in Eq. Q is of the 
form V^(mV^/i), which is irrelevant. In contrast to the 
vapor deposition case, the parameter v can never become 
negative, making a term in unnecessary. Instead, the 
tendency towards phase separation is locally captured by 
K < 0, which is the main focus of the work by Leonard 
and Desai 0, and of Atzmon et al Q (the latter 
study does not include the height profile). A change in 
the sign of K marks the onset of phase separation in 
models of conserved dynamics. Higher order terms, such 
V^m or V^TO^ which are included by other authors, are 
then needed for the stability of short wavelength fluctu- 
ations, and may also affect the precise shape of the order 
parameter profile within domains. Such terms are irrele- 
vant to considerations of the long wavelength behavior of 
Eq. Q), since in this case instabilities can only persists up 
to a length scale of the order ^/WJv, set by the current 
of incoming particles. Since higher temperatures and de- 
position rates favor mixing, it is likely that K eventually 
becomes positive as these parameters are increased. If 
the substrate dimension is d = 1 (or is effectively d — 1 
because diffusion proceeds along a preferred direction) 
the coarse-grained value of K must be positive, as the dy- 
namics are then similar to a 1-dimensional Ising model, 
which cannot have an ordered phase. For this reason, 
the choice of if < in ^ does not capture the long 
wavelength behavior of the system. 

The main focus of the next two sections is on the model 
for vapor deposition, Eqs. Q and (jSJ, and we discuss 
Eqs. (O and only briefly in connection with the RG 
calculation. Our analysis of the models will concentrate 
on the scaling behavior of the height profile and of the or- 
der parameter. On sufficiently large length scales, height 
profiles of growing interfaces are usually characterized by 
a scaling form 



[/i(x,i) - h{x',t')] ) - |x-x' 



'|2X 



\t-t'\ 



(5) 



where x is the roughness exponent, and Zh is a dynamical 
scaling exponent. The values of the exponents depend on 
the underlying growth model, and one of our objectives 
is to find out how they are affected by the coupling to 
the order parameter dynamics. 

The scaling of the order parameter is different along 
the growth direction and perpendicular to it. In contrast 
to the height variable, the order parameter is unlikely 
to be exactly at a fixed point, and for this reason we 
include a correlation length ^. We also have to allow for 
the possibility that the height and the order parameter 
dynamics have different dynamical critical exponents Zfi 
and Zm- The scaling laws for the order parameter then 
read 



Gi:'(x-x') 



(m(x, t)m{'K' ,t)) 



3™(|x-x'|/0 



\t-t'\(^-'y^-gli\t-t'\/C-). (6) 



III. RENORMALIZATION GROUP ANALYSIS 

Let us now renormalize the equations of motion, 
Eqs. and and search for fixed points that are ac- 
cessible by perturbation theory. Inserting the equations 
of motion in the Gaussian probability distribution of the 
noise 

C,,(x,i)' , Cm(x,i) = 



W[Ch,Cm] oc exp |- y d'^x dt 



4D„ 



^ , (7) 

and introducing auxiliary fields m and h, we ob- 
tain the weight of a given space-time configuration 
[ft.(x,t),m(x,i)] 13] 



W[h, m\(x J V[ih] J V[im] exp \j[h, 
with the dynamical functional 



h, fh, m] 



J [ h,h,rh,m]^ J d'^x j dtl^Dh 



hh — h X 



— - lyW^h (Why 



+Dm'mrh — rh x 



dm 
'dt 



KCV^m — rm — uw?) 



~a\/h ■ Vm - hvSj'^h - -m{\/h)'^ 



(8) 



(9) 



The dynamical functional plays the same role in dy- 
namical RG as the Hamiltonian in statics. The bare 
propagators of this model are 



G'Sik,t) = {h{~k,t)hik,t))o = 0{t)e- 



(10) 



Cl}{k,t)^{h{~Kt)hik,t))o^^^^ , (11) 



G^(k, t) = (m(-k, t)m(k, t))o = 0(i)e-^('^+'=')* , (12) 



Co™(k,t) = (m(-k,t)m(k,t))o = 



r) -K{r+k^)\t\ 

K{r + k^) 



(m(x, t)m(x., t')) 



(13) 

and the interaction vertices are obtained from the higher- 
order terms in J. In the diagrams below, h and h are rep- 
resented by straight and wiggly lines respectively. Lines 
for the order parameter m are represented the same way, 
with an additional short dash perpendicular to the prop- 
agator. 
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FIG. 1: The diagram renormalizing Dh. The hnes with a 
small bar represent the order parameter propagators, while 
lines without this bar belong to the height variable. Wiggled 
lines stand for h and m, and smooth lines for h and m. 



The bare dimensions cLq of the couplings are obtained 
by rescaling space, time, height, and order parameter ac- 
cording to X = bx' , t = b^t', h = b^h', h = b^h', m = 
b'^m' , fh = b^fh\ and by requiring invariance of the Gaus- 
sian part of under such rescaUng. (The scahng dimen- 
sion C of the order parameter is related to the exponent 
77 defined in Eq. © via C, = {q — l)/2.) The results are 
listed in table HJ In the following, we analyze the scaling 
behavior resulting from an RG analysis as function of the 
spatial dimension d. 



A. Dimensions d > 6 

In sufficiently high dimensions, the Gaussian fixed 
point, which is characterized by uncoupled, linear 
Langevin equations is stable with respect to the higher- 
order terms. The condition of scale invariance of the 
linear Langevin equations leads to z = 2 and X ^ C = 
(2 — (i)/2, and to the scaling dimensions dg listed in the 
third column of table ^ In dimensions d > 6, all non- 
linear couplings are irrelevant. The surface is smooth, 
and the order parameter goes through a classical phase 
transition. 



B. Dimensions 4 < d < 6 

Below d — the coupling a becomes relevant. The 
Gaussian fixed point still exists, but becomes unstable. A 
new stable fixed point with a nonzero value of a emerges. 
Whenever nonlinear terms cannot be neglected, the cou- 
plings change under rescaling not only according to their 
bare dimensions, but also according to those contribu- 
tions that are generated under renormalization. Renor- 
malization of this model is done by first integrating over 
the large wave vectors K/b < k < A, where A is the 
wave vector cutoff, and the scaling factor b is larger than 
1. Next, the system is rescaled to the original size by 
introducing new variables k' = bk, t' = t/b^. This pro- 
cedure involves an expansion of e"^ in the couplings. In 
this way, the coupling a generates a contribution to Dh, 
which is graphically represented by the diagram in Fig.^ 
Evaluation of this diagram gives a contribution to Dh of 
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do 


dg 


dc 
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b 
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2 - d - 2x 





d - 6 


V 


2-2 








A 


^-2 + X 


(2-d)/2 


4-d 


a 


2C - X + z 


(6-d)/2 







2C - X + ^ - 2 


(2-d)/2 





TABLE I: Bare dimensions do, scaling dimensions dg at the 
Gaussian fixed point, and scaling dimensions da between 4 
and 6 dimensions of all the couplings. 



B = 



A/f)<|k|<A 

a^DlK^A''-%l 



d'^k / dt 





.e~d\ 



{r + A:2)2 



AK^{d-6) 

where Kd is the surface of the d-dimensional unit sphere, 
divided by {^tt)'^. We have also set r = 0, assuming that 
the order parameter is exactly at its critical point. The 
renormalizcd value of this parameter is thus 

D'h = b'-''-^'^[Dh + B]. 

Setting b = I + dl, we obtain the fiow equation 



"^^^ = -Dh{d + 2x-z)+ " 



dl 



4K^ 



(14) 



The exponents z and C are fixed at the values z — 
2 and ^ = (2 — d)/2, since the renormalization of the 
parameters v, K, Dm does not obtain any anomalous 
contributions from diagrams. The condition that a has a 
nonzero fixed point leads to x = 4 — d. With these values 
of the exponents, the condition that Dh is invariant under 
rescaling leads to the fixed point value of Dh 



Dh - 



4/^3(6 - d) 



(15) 



This fixed point, where a is the only nonzero coupling is 
stable between 4 and 6 dimensions. The scaling dimen- 
sions of the other couplings are given in the right-hand 
column of tabled 

Note that for 4 < d < 6, the term arn^ can be regarded 
as a correlated noise acting on the surface height. This 
correlated noise is more relevant than the white noise, 
and incrases the value of the roughness exponent x from 
(2 — c?)/2 to 4 — d. As this value is still negative, the 
surface is flat at this fixed point. 
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(a) /\ k 



(b) 



t=0 





FIG. 2: The diagrams to be considered near d — 4 dimensions. 



C. Dimensions d ~ 4 

Below d — 4 dimensions, the flat phase becomes un- 
stable, because the roughness exponent becomes positive, 
and the coupling A obtains a positive scaling dimension. 
In d = 4-|-e dimensions (with |e| small), we can therefore 
expect a fixed point where A (or a power of A) is of the 
order of e. In order to find this fixed point, let us first 
assume that there is no feedback from the height to the 
order parameter (a = 6 = c = 0), and then take into 
account all terms of the lowest order in A. 

As shown in Figs. |2Ia) and|2Ib), there are two dia- 
grams that contain one A vertex. Diagram (a) makes a 
contribution to q:/2 equal to 



A = 



OO OO 

I df^k I dti I dt2 



A/b<|k|<A ti 



K{r + k^) 

a'^XDrnKiA'dl 

2Kv{K + v) ' 

which is a correction of order Aa^. Fig. \^[h) is a 
correction to Dh of order . Since it modifies the 
flow equation Eq. \i4\ for Dh only to order e, it need 
not be evaluated. Furthermore, Eq. (O describes for 
a ~ h ~ c = Q the relaxational dynamics of an order 
parameter in the universality class of the Ising model, 
which is known to have a non-trivial stable flxed point 
for D.,nu/K = —K4e/9 + O(e^) below 4 dimensions, and 
with r of the order of e 14] . This means that we have to 
take into account additionally diagram|2c) , which makes 
for c? < 4 a contribution 



C = -ZauK J d^k J dt- 

A/6<|k|<A 



Dr. 



'2K(r+k^}t 



eadl 



K{r + B) 



to a/2, which is of order e. In evaluating this expression, 
we have inserted the above-mentioned fixed point value of 
Dmu/K and have set r = 0, considering only the leading 
contribution in an expansion in e. 

Taking all these results together, we obtain the follow- 
ing set of flow equations to order e: 



dl 
dK 

IT 
du 

m 

dX 

n 

dPh 

dl 

da 
H 



Dm{z-d- 20; 
Kiz-2); 
v{z - 2) ; 

\{z-2 + x); 
-Dh{d + 2x- z) 

e 

a 2C~x + z + - 



Kv{K + v) 



(16) 



From the flow equations for K, v and we obtain 
again the fixed point condition z = 2 and C, — (2 — d)/2. 
For e > 0, the flxed point A = is stable, and we have a 
negative roughness exponent x = 4— d, as before. For e < 
0, the fixed point A = is unstable, with the roughness 
exponent x modified due to diagram (c) in Fig. [3 For 
A = and e < 0, Eq. lfTB|l reduces to 



da 

m 



X- 



-X 



2e 



leading to x = 2(4 - d)/3. 

Let us next discuss the fixed point with A 7^ 0. A 
nonzero A requires x — ^- The combination a A then acts 
as an effective coupling, and Eq. H16(l has a non-trivial 
fixed point at 



aX 



Ky{K + v) 

DmK4 



for e > and a fixed point 



2Kiy(K + i^) ,^,0, 
Xa = e ^ — + O(e^), 



(17) 



(18) 



for e < 0. The couplings a, 6, and c have scaling dimen- 
sion zero (to order e) and are thus marginal at this fixed 
point. Determination of their marginal relevance or ir- 
relevance requires evaluation of higher order terms in e, 
which was not attempted in this paper. For e < 0, the 
fixed point is stable as indicated by the flows sketched in 
Fig. 121 We expect it to correspond to a rough phase, with 
the roughness exponent % = 0, possibly receiving correc- 
tions in higher order in e. For e > 0, the fixed point can 
be interpreted as describing a roughening transition. The 
fixed point is unstable (see Figure|31), with flows to either 
a flat phase (if the initial A value is smaller than the fixed 
point value) or to a rough phase, which is not accessible 
by perturbation theory. Compared to a system that is de- 
scribed by the height variable alone, without a coupling 
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FIG. 4: The "brick wall" model used in the simulations. 



FIG. 3: The flow of the coupling Aq near d = 4 dimensions. 
The arrows indicate the direction of the RG flow of Xa. The 
dashed line marks unstable fixed points, and the solid line 
stable fixed points. 



to critical order parameter fluctuations '9] , the lower crit- 
ical dimension for the roughening transition is increased, 
at criticality of order parameter fluctuations, from 2 to 
4. We thus have found that the coupling to the critical 
order parameter fluctuations changes the scaling behav- 
ior of the height variable below 6 dimensions, and that 
it increases the lower critical dimension for the roughen- 
ing transition from 2 to 4. Studying the influence of the 
surface roughness on the critical order parameter fluctu- 
ations near 4 dimensions was not possible to us within 
perturbation theory. To order e, the parameters a, b, 
c are marginal, and to higher orders in e the number 
of diagrams becomes large. Furthermore, it is not clear 
whether further fixed points except those discussed by us 
are at all accessible by perturbation theory. 

In the case of the conserved Langevin Eqs. Q and Q 
proposed for MBE growth, the mutual effects between 
height and order parameter are much weaker. The Gaus- 
sian fixed point is stable above d = 2 dimensions. Below 
d — 2, the coupling /3 becomes relevant. It K > and v 
is small {v <C A'^/K), the diffusion of the order parame- 
ter on the surface affects the height profile. Performing 
an analysis analogous to the one above near 6 dimensions 
(where a was the relevant coupling), we find a stable fixed 
point in d = 2 - e with f3^ = DhK^TTe/2Dl^ + 0{e^) and 
X = 2 — d. (Without coupling to the order parameter, 
there is a smaller roughness exponent of x = (2 — d)/2.) 
For negative K, there can be other interesting (non crit- 
ical) effects, as mentioned before. 



IV. COMPUTER SIMULATIONS 

While the RG gives information about scaling behavior 
in high dimensions (4 and higher), computer simulations 
are particularly efficient in low dimensions. In this sec- 
tion, we present results from simulations in 1+1 dimen- 
sion. Rather than discretizing Eqs. 11I2|) and integrat- 
ing them numerically, we perform numerical studies of a 
"brick wall" restricted solid-on-sohd model (see Fig. 
Since this model shares the same symmetries and con- 
servation laws as the Langevin equations, it should share 
some of the same universality properties. However, as 
noted by Oskoee, Khajehpour, and Sahimi there are 



also important distinctions between the continuum and 
discrete model. In the continuum model at zero tempera- 
ture, the coarsening force on a domain of length £ decays 
exponentially with £. The corresponding coarsening time 
then grows logarithmically in £ (z oo) as opposed 
to z = 2 for Glauber dynamics [l3 |. 

Starting from a flat surface, particles are added such 
that no overhangs are formed, and with the center of each 
particle atop the edge of two particles in the layer below. 
We use two types of particles, A and B (black and grey 
in the figure). The probability for adding a particle to 
a given surface site, and the rule for choosing its color, 
depend on the local neighborhood. Since growth is slower 
on slopes, these growth rules correspond to A < [T9ll20l |. 

When A (B) particles are more likely to be added to A 
dominated regions (B), the particles tend to phase sepa- 
rate and form domains. In this case, the order parameter 
correlation length ^ is of the order of the average domain 
width. By changing the growth rules, it is possible to 
study cases in which some (or all) of the couplings a, b, 
c, and a vanish, and thus to gain a more complete picture 
of the different ways in which the height and the order 
parameter influence each other. 

When all the couplings between the order parameter 
and the height vanish {a = b = c — a = 0), the well- 
known critical exponents Zh = 3/2 and x = 1/2 of the 
KPZ equation |9Ml5j, and z„i = 2 and = 1 of the 
Glauber model 'l3| are recovered. This situation is im- 
plemented in the following way: A surface site is chosen 
at random, and a particle is added if it does not gener- 
ate overhangs. Its color is then chosen depending on the 
colors of its two neighbors in the layer below. If both 
neighbors have the same color, the newly added parti- 
cle takes this color with probability 1 — p, and the other 
color with probability p (where p is much smaller than 1). 
If the two neighbors have different colors, the new par- 
ticle takes either color with probability 1/2. Neighbors 
within the same layer are not considered. As discussed in 
Ref. 12], these growth rules lead to an order parameter 
correlation length ^ ^ l/-\/P as p — s- 0. 

Here, we want to focus on the more interesting situa- 
tions where either a or a, 6, c (or all of them) are nonzero. 
In the first case presented below, the order parameter af- 
fects the height variable, but is not influenced by it. In 
the second case, the height profile affects the order pa- 
rameter dynamics, but not vice versa. One would expect 
that in the first case the order parameter imposes its dy- 
namical exponent z = 2 on the dynamics of the height 
profile, and that in the second case the height profile im- 
poses its exponent z = 3/2 on the order parameter. This 
latter is, however, not the case for a/A < 0, and we shall 
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FIG. 5: The scaling function g{y), obtained for L = 2048 by 
collapsing data for p = 320, 640, 1280, and 2560. For r, the 
values 0.05 and 0.025 were used. The dashed line is a power 
law ocl/y. 

see that Zm is nonuniversal in this case. In the fuUy cou- 
pled case with a, a, b, and c not equal to zero, we find 
z = 2orz = 3/2, depending on the sign of (Aa). Some 
of these results were already reported in Ref. 

A. Growth influenced by independent phase 
ordering (a 7^ 0, a, 6, c = 0) 

The situation a > (a < 0) is implemented by updat- 
ing sites on top of particles of different colors less (more) 
often by a factor r < 1 (r > 1) compared to sites above 
particles of the same color. If the color of the new particle 
depends only on the neighbors in the layer below, the or- 
der parameter is not affected by the height variable, and 
its dynamics is still the same as that of an Ising model, 
with = 2. 

We first discuss the case a > 0: Because growth is 
slower at domain boundaries than within domains, the 
domain boundaries sit preferentially at the local minima 
of the height profile, with a mound over each domain. 
This implies that the surface roughness exponent is x = 1 
on length scales up to ^. Changes in the height profile 
on this scale result from domain wall diffusion, and the 
dynamical exponent is therefore — 2. On length scales 
much larger than ^, the average order parameter is zero, 
implying that KPZ exponents of x = 1/2 and Zh = 3/2 
are regained. The crossover in the roughness is described 
by a scaling form 

{[h{x, t) - h{x\ m = \x- x'l'g (^^-^^ , 

with a constant g{y) for ?; ^ 1, and g{y) ~ 1/y for 
y ^ 1. Figure El shows our simulation results for g{y), 
obtained by the data collapse of {[h{x, t) — h{x', t)]^)/|x — 
x'l'^ versus \x—x'\y/p, for different values of p. The curves 
for the two different values of r differed only slightly, and 



were collapsed by multiplying the curves for r = 0.05 by 
a factor of 1.06. The scaling collapse is compatible with 
g{y) ^ 1/y for large y, and g{y) — s-constant, for small y. 

For a < 0, growth occurs with a larger probability at 
domain boundaries. Therefore, domain boundaries sit at 
local maxima. However, further away from the domain 
boundaries, their effect is not felt, and we find % = 1/2 
and Zh = 3/2, just as in the case a = 0. 



B. Phase ordering influenced by independent 
growth (a = 0; a, 6, c 7^ 0) 

The situation a = is implemented by choosing r — 1, 
i.e., adding a particle at each possible site with the same 
probability, irrespective of the color of its neighbors. To 
mimic the influence of surface roughness on the order 
parameter (nonzero a, b, or c in Eqs.Q), the color of 
a newly added particle is made dependent not only on 
those of its two neighbors in the layer below, but also on 
the colors of its two nearest neighbors on the same layer, 
if these sites are already occupied. With probability 1—p, 
the newly added particle takes the color of the majority of 
its 2, 3, or 4 neighbors, and with probability p it assumes 
the opposite color. If there is a tie, the color is chosen 
at random with equal probability. Since the neighbor on 
the hillside of a site is more likely to be occupied than 
the one on the valley side, with this rule domain walls 
are driven downhill, corresponding to a > in Eq. 1^)1. 
Also, domains on hilltops can expand more easily than 
those on slopes or in valleys, suggesting a value of 6 > 0. 

We reported already in Ref. [iJI that the dynamical 
critical exponent Zm associated with the order parameter 
has a nontrivial dynamic exponent of z^ — 1.85, and not 
the value z^ = 1-5, which may be expected if the walls 
follow the surface fluctuations. A potential explanation 
was provided in Ref. in connection with the dynam- 
ics of a single domain wall riding on a growing surface. 
The growth rules for the surface imply that sequences of 
brick addition usually proceed from local minima in the 
uphill direction, since the addition of a brick generates a 
potential growth site (where a brick can be added with- 
out generating overhangs) at the nearest uphill position. 
The walls that try to slide downhill are therefore faced 
with an upward avalanche of growth mounds of different 
sizes that hamper their downhill motion. 

The exponent Zm is not only nontrivial, but also 
nonuniversal and depends on the value of a. A change 
in a can be implemented in the computer simulations by 
taking into account neighbors within the same layer as 
the site being updated with a probability q smaller than 
1. Using the above-mentioned rules and values of q = 1, 
0.25, 0.125, 0, we find z„ = 1.8, 1.89, 1.96, 2.0. The 
difference from the value for g = 1 reported in Ref. 
stems from the fact that in that paper we had inadver- 
tently assigned to neighbors within the same layer a dou- 
ble weight - thus illustrating once more the nonuniversal- 
ity of the value of Zm- Each of these values was evaluated 
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FIG. 6: Scaling collapse of the correlation functions Gm for 
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FIG. 7: Domain density (i.e., number of domains, divided 
by the system size) as a function of time for L — 16384, 
averaged over 100 samples, for 4 different values of q. The 
curves are well fitted by power laws of the form t^^/^™ with 
the exponents Zm=l-S, 1.89, 1.96, 2.0. 



in two independent ways, in order to confirm its stabil- 
ity and parameter dependence. The first was a collapse 
of the correlation function Gm^ (x) , as shown in Fig. (HJ 
Since p sets the inverse time scale, the scaling variable 
is yielding the value of Zm- The second method 

was a domain coarsening simulation, following a quench 
from p — 0.5 to p = 0. Fig. |7| shows the number of do- 
mains, divided by the system size, as function of time for 
different values of q. 

Equations H1I2|) allow also for the possibility of a < 0. 
This would imply that domain walls move uphill and 
that identical neighbors in the same layer have a positive 
interaction energy, in contrast to neighbors in different 
layers. Although this is an implausible physical situa- 
tion, it is nevertheless of some theoretical interest. For 
a = X and c = 0, the invariance under the transforma- 
tion x' = X — Xet and h' — h + ex (with a small e) of 
the KPZ equation holds also for Eq. This invariance 



corresponds to the Galilean invariance of Burger's equa- 
tion. It implies that the order parameter dynamics are 
goverened by the same time scale as the height variable, 
i.e., Zm = Zh- Implementing the case a < in our sim- 
ulations, we indeed find z,„ = 1.5, suggesting that the 
Galilean invariant fixed point is the only attractive fixed 
point in this domain of parameters. 



C. Mutual couplings (a, a, b, c ^ 0) 

When all couplings are different from zero, the proba- 
bility for adding a particle at a given site and the choice 
of the particle color depend on the local neighborhood. 
The simulation parameter q is positive, and r is different 
from 1. 

For a > 0, we find Zm = '2 irrespective of the values of 
a, b, and c. As particles are added to domain boundaries 
with a smaller probability than within domains, domain 
boundaries sit at local minima most of the time. There- 
fore they perform a random walk even when a ^ 0. Over 
each domain there is a mound, implying that x — 1 on 
length scales up to ^. Changes in the height profile on 
this scale result from domain wall diffusion, and the dy- 
namic exponent is therefore z/j = 2. On length scales 
much larger than ^, the average order parameter is zero, 
and KPZ exponents of x = 1/2 and Zh — 3/2 are re- 
gained. 

For a < 0, particles are added rapidly on domain 
boundaries, and domain walls can therefore not be 
trapped in local height minima. For the case aX > 0, 
where domain walls tend to move uphill, we therefore ex- 
pect that the situation a < is similar to the case a = 0, 
for which we found z^, = Zh = 5/2. This means that the 
above-mentioned Galilean-invariant fixed point remains 
applicable to a < 0. The simulation results shown in 
Fig. |S1 confirm this expectation. For the case aX < and 
a = 0, we have argued above that the downhill motion 
of domain walls is hampered by an upward avalanche of 
growth mounds, which cause the walls to be temporarily 
stuck in local minima, leading to a nonuniversal exponent 
Zm- Now, for a < and aX < 0, we find in our simula- 
tions that the dynamical critical exponent z,„ is identical 
to Zfi = 1.5, implying that the downhill motion of the 
domain walls is not hampered any more but that the 
domain walls can follow the height fluctuations. Fig. |S1 
shows the results of a domain coarsening simulation for 
the parameters q = 1 and r = 5. For comparison, sim- 
ulation results for positive aX (and otherwise the same 
parameter values) are also shown. One can see that the 
exponent Zm is indeed the same in both situations. To 
summarize, we find that for q < the height profile 
imposes its critical behavior on the order parameter fluc- 
tuations, while in the opposite case, a > 0, the domain 
wall diffusion imposes a dynamical exponent z = 2 on 
the system. 
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FIG. 8: Domain density (i.e., number of domains, divided 
by the system size) as a function of time for aX < and 
L = 16384, averaged over 100 samples, for q = 1 and r — 5 
(solid line). The dashed line is a power law with an exponent 
—2/3, corresponding to Zm = 1-5. The dotted curve shows 
the corresponding simulation result with positive aX. 



D. Comparison to analytical results 

It is interesting to compare the results of the computer 
simulations with the admittedly limited results of the RG 
analysis presented in the previous section. 

Let us start from our last result that for aA > the 
height profile imposes its critical behavior on the order 
parameter, while the reverse is true for aX < 0. The RG 
up to order e = 4 — d showed that for aA < the fixed 
point is accessible perturbatively, with z = 2 to order e. 
It is striking that this result holds also in 1 + 1 dimension. 
For q;A > 0, the RG flow runs away to infinity, suggesting 
the existence of a strong coupling fixed point. The critical 
behavior corresponding to this fixed point was found in 
our simulations in 1+1 dimension to be the same as that 
of the KPZ equation, with z = 3/2. If this result is not 
particluar to 1+1 dimension, it suggests that a positive 
aA > might be irrelevant at the KPZ strong coupling 
fixed point. However, this conclusion can not be tested 
via by perturbative RG analysis. 

Up to order e of the RG analysis, the parameters a, 6, c 
are marginal. It appears from our simulations that these 
three coupling indeed do not modify the critical behavior 
as long as a 7^ 0, suggesting that these parameters are 
marginally irrelevant. 

For a ~ 0, the critical behavior of the height profile is 
given by the KPZ equation, and it has a stable fixed point 
at A = in which the dimension of ft, is (2 — d)/2 < 0. 
Thus at the weak coupling fixed point, a, b, and c are 
irrelevant in four dimensions and the Ising fixed point is 
not modified by coupling to the height parameter. How- 
ever, the RG analysis cannot predict the critical behavior 
of the coupled system in the rough phase, which is char- 
acterized by a strong coupling fixed point. Nevertheless, 
we could argue that there exists a Galilean-invariant fixed 
point when aA > 0, where the height profile imposes its 



dynamical critical exponent on the order parameter. The 
result Zm = 2 in 1+1 dimensions was confirmed by our 
computer simulations, suggesting a larger domain of pa- 
rameter space where such scalings apply. 

The only case for which we have no analytical predic- 
tion is when a — 0, and Xa < 0, where the computer 
simulations reveal nonuniversal behavior. 



V. RELATION TO PASSIVE SCALAR 
ADVECTION AND DRIFTING POLYMERS 

There is a close connection between the model stud- 
ied in this paper and other coupled nonequilibrium sys- 
tems. One such system is obtained when we regard do- 
main walls as "particles" that ride on the growing sur- 
face. With the substitution Vto — p, we obtain for 
u = b = c = the following equation for p: 

1^ =AcV2p + aV(pV/i) + Cp(x,t), (19) 

with a conserved noise ^p. Without the noise term this 
equation is the Fokker Planck equation corresponding to 
the Langevin equation 

H-x 

^aVh + Ut). (20) 

at 

If we assume that the particles do not interact with each 
other, the Fokker Planck equation, combined with a noise 
term, describes the time evolution of the particle density. 

If we require additionally that there is no effect of the 
particles on the growing interface, the coupling a van- 
ishes, too. The substition V/i = v then turns the KPZ 
equation into a randomly stirred Burger's equation, and 
Eq. (|20|l becomes the equation of motion of a particle ad- 
vected by the flow. This model was studied in detail in 
Ref. |2l|. Just as for the model described in this paper, 
the scaling behavior of the advected particles is funda- 
mentally different in the two cases a/A > and a/ A < 0. 
In the first case, the system has a Galilean invariant fixed 
point, and particle diffusion is characterized by a dynam- 
ical critical exponent Zp = 3/2 in one dimension, while 
this exponent is larger than 3/2 and nonuniversal in the 
other case. 

Finally, Eqs. I|1I2|) with a = b and c — r = u = 0, 
but with a ^ 0, can be mapped on the equations used 
to describe the dynamics of a streched string moving in 
a random medium |2^ . If the string is streched in the 
x-direction, and if it moves in the /i-direction, the con- 
figuration of a string embedded in 3 dimensions can be 
characterized by the coordinates h{x) and h±{x). As- 
suming that the evolution of the line is dissipative and 
local, the equations of motion then are our equations 
Eqs. (|1I2|I with the replacement m — dxh± and with Cm 
replaced with a conserved noise dxCm- Apart from the 
Galilean- invariant fixed point, this set of equations has 
a fixed point where a fluctuation-dissipation relation is 
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satisfied (i.e. where a stationary solution of the Fokker- 
Planck-equation can be written down). This fixed point 
corresponds in our notation to the situation a = Ka/v. 
(For a general discussion of equations the stationary so- 
lutions of which can be calculated, see Ref. 0). 

VI. CONCLUSIONS 

In summary, the interplay between surface roughening 
and phase separation leads to a variety of novel criti- 
cal scaling behaviors. At one extreme, the height profile 
adapts to the dynamics of critical domain ordering. At 
the other extreme, the dynamics of the domain walls fol- 
low the height fluctuations. For a third range of paramter 
values, the dynamics of domain wall motion is influenced 
by the roughness, but exhibits nontrivial and nonuniver- 
sal scaling behaviors. 

Several generalizations of the model presented in this 
paper are possible. For example, as discussed in Ref. p3| . 
one can consider the situation where the symmetry 
breaking involves a continuous, rather than an Ising-like 
order parameter. Such a situation applies to the depo- 
sition of spins that can realign on the surface but are 



frozen in the bulk, or to orientational symmetry break- 
ing in the plane during crystal growth. Another gener- 
alization would be the inclusion of elastic forces, which 
are often present during the growth of composite films 
(see Ref. |2J|). Furthermore, one could consider phase 
transitions where the different types of molecules order 
on sublattices instead of phase separating. 

Finally, there exist growth situations where one type of 
particles is magnetic. In addition to a ordering or phase 
separation transition which occurs at the surfaces, there 
is in this case also a magnetic phase transition, which 
occurs in the bulk. This combination of two-dimensional 
and three-dimensional phase transitions is particularly 
challenging for a theoretical analysis, and it leads to in- 
teresting experimental results p5| . 
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